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In this Letter we report Ising model simulations of the growth of alloys which predict quite different 
behavior near and far from equilibrium. Our simulations reproduce the phenomenon which has been 
termed “solute trapping,” where concentrations of solute, which are far in excess of the equilibrium 
concentrations, are observed in the crystal after rapid crystallization. This phenomenon plays an 
important role in many processes which involve first order phase changes which take place under 
conditions far from equilibrium. The underlying physical basis for it has not been understood, but these 
Monte Carlo simulations provide a powerful means for investigating it. 


PACS numbers: 61.50.Cj, 61.46.+W, 64.60.-i 


Rapid cooling is being used increasingly to produce 
fine, tailored microstructures, to produce fine dispersions 
by rapid quenching of small particles, and to modify the 
surface properties and structure of materials by rapid ther- 
mal processing or by laser irradiation. In some materials, 
the amorphous structure of the liquid is preserved in the 
form of a metastable glassy phase [1-3]. In other ma- 
terials, equilibrium phases are suppressed, and unstable 
or metastable phases, which do not exist otherwise, form 
during a rapid quench. In most materials, compositions 
outside the equilibrium range can be obtained. Related 
nonequilibrium segregation effects such as the “facet ef- 
fect” in the Czochralski growth of silicon [4,5] are present 
at relatively slow growth rates. There is no sound phys- 
ical understanding or theoretical framework for these and 
many other nonequilibrium effects associated with first or- 
der phase changes which take place far from equilibrium. 

The best data available for the formation of phases 
which have nonequilibrium compositions have been ob- 
tained for the incorporation of dopants into silicon during 
very rapid solidification following laser melting of a sur- 
face layer [6-14]. These dopants are incorporated into 
the crystal at high concentrations which are metastable at 
any temperature, and they will precipitate during subse- 
quent annealing. The quantitative measure of this effect 
is the distribution coefficient, also known as the k value, 
which is the ratio of the concentration of dopant in the 
solid at the interface to the concentration of dopant in 
the liquid at the interface. The k value increases dra- 
matically with growth rate in these experiments, from 
an equilibrium value of 0.001 to a value approaching 1. 
Several models have been proposed to account for the 
phenomenon [15-18], and the model which has been 
compared most extensively with experiment is due to 
Aziz [17-21]. The underlying physical basis for this 
model is controversial, but it does provide a distribution 
coefficient which can differ significantly from the equi- 
librium distribution coefficient, and it has been used to 


fit to experimental results using reasonable values for the 
diffusion coefficient. 

Previous Monte Carlo simulations, based on the Ising 
model applied to alloys, have not reproduced these ex- 
perimental results [22]. Why this is so has been a mys- 
tery, since our present understanding of the atomic scale 
processes involved in crystallization is based on the Ising 
model, and Monte Carlo simulations [23-26] guided the 
development of this understanding, which includes the sur- 
face roughening transition and its central role in the equi- 
librium properties of surfaces and interfaces, and in the 
kinetics of interface motion during crystal growth [27]. 

The simulations reported here have been done in two 
dimensions and are a special version of a two-dimensional 
random walk. In order to simulate crystal growth in two 
dimensions, the position of the “walker' becomes a line 
which represents an interface. In three dimensions, a 
surface represents the interface. The crystal is on one side 
of this line or surface, and the liquid is on the other. The 
probability that the interface jumps backward or forward 
at a site depends on how many of its nearest neighbors 
are on the crystal side or on the liquid side of the 
interface. This lateral coupling tends to keep the interface 
planar, and so effectively introduces a surface tension. 
The probability of interface jumps in the two directions 
can be biased to simulate nonequilibrium melting or 
freezing. To model an alloy, special sites are introduced 
to represent the alloying element or dopant. There is a 
different bias for interface jumps across these sites than 
at the normal sites. The simplest case conceptually is 
when these special sites are stationary. This corresponds 
to diffusionless growth, which is the limiting case for 
the diffusion of atoms being slow compared to the rate 
at which the interface moves. In this limit, the liquid 
and solid both have the same composition since only the 
interface moves, not the atoms, and so the k value is 
1. The net rate of motion of the interface depends on 
an average of the biases for the two types of sites. The 
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special sites can also be made to move around in order to 
simulate diffusion. In most alloys, diffusion in the liquid 
is very rapid compared to diffusion in the solid, and so it 
is usual to assume that there is diffusion only in the liquid. 
The usual near-equilibrium growth conditions for alloys 
are reproduced when the special sites move around rapidly 
compared to the rate at which the interface moves. The 
“equilibrium" k value in this case depends on the relative 
biasing for the two types of sites. When the biasing is 
adjusted so that the equilibrium k value is less than 1, the 
dopant atoms tend to stay in the liquid and diffuse away 
from the interface. But, as pointed out above, the k value 
is 1 if the special sites are stationary, independent of the 
biasing. There is a transition between these two extremes 
which occurs when the jump rate for the special sites is 
comparable to the rate at which the interface moves. We 
believe this to be the origin of solute trapping behavior, 
and it is evident that this effect should be present in this 
model as well as in experiment. 

Simulations were reported earlier [28,29] using this 
model for “diffusionless” transformations for which the k 
value is 1. These provided the first clear-cut confirmation 
of the theoretically expected behavior for a diffusionless 
transformation, that is, that “freezing” or “melting” should 
be reversible for temperatures above and below T 0 , the 
(composition dependent) temperature at which the free 
energy of a solid alloy is equal to the free energy of 
a liquid alloy with the same composition. This is a 
significant result because the jump rates of the interface 
were specified in terms of the chemical potentials of 
the atoms, not in terms of the free energies of the 
phases. This suggests that the simulations contain the 
correct physical model to explain solute trapping. These 
simulations also provided the first information about 
crystallization kinetics in diffusionless transformation. 

The Monte Carlo simulations reported here are in the 
transition region where the k value is intermediate be- 
tween the equilibrium value and 1, and the diffusion coef- 
ficient and the growth temperature in the simulations have 
been adjusted to reproduce the experimentally measured 
distribution of a dopant in a silicon crystal after rapid re- 
crystallization of a surface layer. The transition rates at 
the interface for the Monte Carlo scheme used here are 
similar to those used previously [22,28,29]. The transi- 
tion probability for an atom to leave the crystal to join the 


liquid can be written as 

Pk = p °k exp 


V[ $ki_ 

^ k B T 


( 1 ) 


Here P° is a constant, k B is Boltzmann’s constant, and T 
is the temperature. The superscript C identifies the atom 
as being in the crystal, and the subscript k identifies the 
species of the atom. The summation is over the nearest 
neighbors of the atom; each nearest neighbor is identified 
by a pair of indices j and /, where j is either C (crystal) or 
L (liquid), and / identifies the species of the neighboring 


atom. is the energy of the bond between one atom 
defined by ik and its neighbor defined by jl. Similarly, 
the probability of an atom going from liquid to crystal is 


p: = p°l 


k exp 



( 2 ) 


where AS* is the entropy of fusion for the species k. This 
form preserves microscopic reversibility. The simulations 
reported here start with both solid atoms and liquid atoms 
present in a two-dimensional array with an interface 
between them. Diffusion in the liquid is modeled by the 
interchange of two adjacent atoms of different species, 
with a probability which is chosen so that the average 
jump rate T is some multiple (or fraction) of P*. A 
temperature is chosen, and the transition probabilities P* 
are calculated from Eqs. (1) and (2). Individual atoms 
at the interface join or leave the crystal based on these 
transition probabilities. A normalized rate of motion of 
the interface is then calculated from the net motion of the 
interface. 

Figure 1 shows typical experimental data for silicon 
implanted with bismuth [10]. The as-implanted distribu- 
tion of bismuth is shown as the dashed line. The sam- 
ple was surface melted with a laser to a depth of about 
2 /im. The interface layer stayed molten for about 2 /zs, 
during which time the implanted atoms diffused in the 
liquid, spreading out the as-implanted distribution. The 
interface then came back toward the surface, pushing the 
bismuth atoms ahead of it, resulting in the distribution de- 
termined by Rutherford backscattering as shown by the 
open circles. An experimental distribution coefficient ( k 
value) was determined [10] by estimating the amount by 



FIG. 1. The open circles are the experimental data for the 
depth distribution of bismuth implanted into silicon after 
laser melting of a surface layer. The closed circles are the 
depth distribution of bismuth atoms from the Monte Carlo 
simulations. 
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which the implanted bismuth distribution spread out by 
diffusion while the surface layer was liquid, then adjust- 
ing the k value to fit the final bismuth distribution which 
resulted from the passage of the interface. Although the 
equilibrium k value for bismuth in silicon is 8 X 10~ 4 , 
the data were fitted with a k value of 0.1. This ex- 
periment, and many other similar experiments, unequiv- 
ocally established that the k value is strongly growth rate 
dependent. 

Simulations were performed for dilute Si:Bi alloys by 
using bond energies, which were calculated from the 
thermodynamic properties of the alloy using a regular 

7g AgsiAg ASeiAa k t 
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solution model for the solid and an ideal solution model 
for the liquid, with an equilibrium k value of 8 X 10 -4 . 
The entropy of fusion, which determines the roughness of 
the interface at equilibrium, was chosen to stay below the 
two-dimensional critical point, and to make the repeatable 
step density, which depends on the roughness of the 
interface, similar to that for silicon. A repeatable step site 
(also known as a kink site) is an interface site which has 
half of its nearest neighbor sites occupied by atoms of the 
crystal. The same values were used for the CL bonds and 
for the LL bonds. The values used in these simulations 
are listed here: 

(<£sisi 0SiSi )/ka (^B^i “ \)/ka ($SiBi ~ <AsiBi )As 

8440 K 2722.5 K 5430.8 K 


Tm 10(1 T m are the melting points of silicon and bismuth, 
respectively. 

A two-dimensional array of 150 X 1500 atoms was 
used, which corresponds to a sample that is about 45 nm 
wide by 450 nm deep. This depth is about the same 
as the total penetration depth of bismuth atoms in the 
experiment. Since the as-implanted profile spreads out 
due to diffusion in the liquid during the time that a 
liquid layer exists, a profile for the distribution of bismuth 
atoms at the time when the crystallization front reached 
it was established as follows. 63 bismuth atoms were 
inserted into a layer 140 nm deep in the crystal, which 
corresponds to the total implant concentration and to the 
average implant depth. These atoms were then allowed 
to diffuse to a half- width of the order of 100 nm, 
comparable to the width of the bismuth distribution 
in the experiment just before the recrystallization front 
reached the deepest bismuth atoms. The same distribution 
could have been achieved by letting the atoms diffuse 
while the interface penetrated into the crystal and then 
returned to the depth of the bismuth atoms, but the 
crystallization front does not interact with the bismuth 
atoms during this period. The actual rms half-width 
of the distribution used was 340 atomic layers. The 
closed circles in the figure show the distribution of the 
bismuth atoms in the simulations after the passage of the 
crystallization front at a temperature of 1400 K, with the 
diffusion jump rate 40 times the arrival rate of atoms at a 
repeatable step site. The final positions of the atoms were 
smoothed by using Gaussian averaging to simulate the 
experimental broadening associated with the Rutherford 
backscattering depth measurement. The wiggles in the 
final distribution results from having a total of only 
63 bismuth atoms in the simulation. The simulations 
reproduce the experimental result, giving a k value of 
0.1 for these growth conditions. Similar agreement has 
been obtained for data in Ref [13]. Simulations run at 
other temperatures (growth rates) and with various liquid 
diffusion coefficients indicate that the k value depends on 


both. Similar results have been obtained with simulations 
in three dimensions. 

In the previous Monte Carlo simulations [22], atoms 
arrived randomly at surface sites, as given by Eq. (2), 
and departed from surface sites as given by Eq. (1). The 
k value was taken as the ratio of the concentration of 
dopant in the growing crystal to the concentration of 
dopant in the arriving atoms. The arrival and depar- 
ture rates from the crystal were identical to those used 
here, but there was no liquid phase present. This was be- 
lieved to be a valid model since dopant incorporation into 
a crystal, including solute trapping, is a process which 
takes place at the interface and provides the boundary 
condition at the interface for the diffusion field in the 
liquid. This implies that the diffusion process in the liq- 
uid can be separated from what happens at the interface. 
However, in these simulations [22] it was found that the 
k value did not increase very much from the equilibrium 
value with growth rate. Agreement with experimental 
data could be obtained only by introducing surface seg- 
regation and stress near the interface. The results pre- 
sented here indicate that the essential difference between 
the two simulation schemes is the presence of the liquid 
phase. 

The Ising model Monte Carlo computer simulations 
reported here reproduce experimental results on solute 
trapping with the only input being equilibrium thermo- 
dynamic data, a diffusion coefficient for the atoms in the 
liquid phase, and a growth temperature. This model is 
actually the simplest possible scheme for simulating al- 
loy crystallization, but it exhibits unexpected and quite 
complex behavior. We plan to explore this behavior in 
order to develop an understanding of the implications 
of the simple underlying assumptions. These simula- 
tions provide a powerful and flexible means of exploring 
the consequences of these simple assumptions for crystal 
growth, and it is expected that this will provide significant 
new insights into the crystallization of multicomponent 
materials. 


2532 


Volume 75, Number 13 


PHYSICAL REVIEW LETTERS 


25 September 1995 


The part of this work which was carried out at the 
University of Arizona was supported by NASA Contract 
No. NAG8-944. 


[11 P. Duwez, ASM Trans. Q. 60. 606 (1967). 

[2] P. Duwez. Prog. Solid State Chem. 3. 377 (1967). 

[3] Ultrarapid Quenching of Liquid Alloys, edited by 
H. Herman (Academic Press, New York, 1981). 

[4] T. Abe, J. Cryst. Growth 24/25. 463 (1974). 

[5] A. J. R. de Kock, in Crystal Growth and Materials, edited 
by E. Kaldis (North- Holland, Amsterdam. 1977), p. 693. 

[6] P. Baeri, J. M. Poate, S. U. Campisano. G. Foti. E. Rimini, 
and A.G. Cuilis, Appl. Phys. Lett. 37, 912 (1980). 

[7] P. Baeri, G. Foti, J. M. Poate, S. U. Campisano, and A. G. 
Cuilis, Appl. Phys. Lett. 38, 800 (1981). 

[8] C.W. White, S. R. Wilson, B.R. Appleton, and F. W. 
Young Jr., J. Appl. Phys. 51, 738 (1980). 

[9] C.W. White, B.R. Appleton, B. Stritzker. D M. Zehner. 
and S.R. Wilson, Mater. Res. Soc. Symp. Proc. 1, 59 
(1981). 

[10] P. Baeri, G. Foti. J. M. Poate, S. U. Campisano, E. Rimini, 
and A.G. Cuilis, Mater. Res. Soc. Symp. Proc. 1, 67 
(1981). 

[11] C.W. White, D M. Zehner, S. U. Campisano, and A.G. 
Cuilis, Surface Modification and Alloying (Plenum. New 
York, 1983), p. 81. 

[12] M.J. Aziz, J.Y. Tsao, M.O. Thompson, P. S. Peercy, 
C. W. White, and W. H. Christie, Mater. Res. Soc. Symp. 
Proc. 35, 153 (1985). 

[13] M.J. Aziz, J.Y. Tsao, M.O. Thompson. P.S. Peercy, and 
C.W. White, Phys. Rev. Lett. 56, 2489 (1986). 


[I4| M.J. Aziz and C.W. White. Phys. Rev. Lett. 57. 2675 
(1986). 

[15] K. A. Jackson. G.H. Gilmer, and H.J. Leamy. in Laser 
and Electron Beam Processing of Materials, edited by 
C.W'. White and P.S. Peercy (Academic. New York. 
1980), p. 104. 

[16] R.F. Wood. Appl. Phys. Lett. 37. 302 (1980). 

[17] M.J. Aziz, J. Appl. Phys. 53. 1158 (1982). 

[18] M.J. Aziz, Appl. Phys. Lett. 43, 552 (1983). 

[19] M.J. Aziz. Mater. Res. Soc. Symp. Proc. 23. 369 (1984). 
[201 L.M. Goldman and M.J. Aziz. J. Mater. Res. 2. 524 

(1987). 

[21] M.J. Aziz and T. Kaplan, Acta Metall. 36. 2335 (1988). 

[22] G.H. Gilmer. Mater. Res. Soc. Symp. Proc. 13. 249 
(1983); Mater. Sci. Eng. 65. 15 (1984). 

[23] H.J. Leamy and K. A. Jackson. J. Appl. Phys. 42, 2121 
(1971). 

[24] H.J. Leamy and G.H. Gilmer, J. Cryst. Growth 24/25. 
499 (1974). 

[25] H.J. Leamy, G.H. Gilmer, and K. A. Jackson. Surface 
Physics of Materials, edited by J.M. Blakeley (Academic. 
New York, 1975), p. 319. 

[26] G. H. Gilmer and K. A. Jackson. 1976 Crystal Growth and 
Materials, edited by E. Kaldis and H.J. Scheel (North- 
Holland. Amsterdam, 1977), p. 80. 

[27] J.D. Weeks and G. H. Gilmer, Adv. Chem. Phys. 40. 157 
(1979). 

[28] K.A. Jackson, G.H. .Gilmer, D. E. Temkin. J.D. Wein- 
berg. and K. Beatty, J. Cryst. Growth 128, 127 (1993). 

[29] K. A. Jackson, G.H. Gilmer, D. E. Temkin. and 
K. Beatty, in Proceedings of the East- West Surface 
Science Workshop '94 (to be published). 


2533 



